<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of AM_example</title>
  <meta name="keywords" content="AM_example">
  <meta name="description" content="">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html v1.5 &copy; 2003-2005 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../index.html">Home</a> &gt;  <a href="index.html">CO4OI</a> &gt; AM_example.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../index.html"><img alt="<" border="0" src="../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="index.html">Index for CO4OI&nbsp;<img alt=">" border="0" src="../right.png"></a></td></tr></table>-->

<h1>AM_example
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong></strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre class="comment"></pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
This function is called by:
<ul style="list-style-image:url(../matlabicon.gif)">
</ul>
<!-- crossreference -->



<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../up.png"></a></h2>
<div class="fragment"><pre>0001 clear;
0002 
0003 
0004 addpath(genpath(<span class="string">'misc'</span>))
0005 addpath(genpath(<span class="string">'operators'</span>));
0006 addpath(genpath(<span class="string">'Solvers'</span>));
0007 addpath(<span class="string">'Data'</span>)
0008 
0009 <span class="comment">%% Generate image example</span>
0010 
0011 <span class="comment">% input parameters</span>
0012 n= 16; <span class="comment">% n rows, image dimension n x m</span>
0013 m= 16; <span class="comment">% m columns</span>
0014 N=n*m;
0015 plots=1; <span class="comment">%auxiliary flag for plots</span>
0016 noisy=1; <span class="comment">%auxiliary flag indicating presence of noise. 1: noisy case; 0: noisless case</span>
0017 input_snr=30; <span class="comment">%input snr</span>
0018 
0019 
0020 
0021 <span class="comment">% parameters image</span>
0022 <span class="comment">% field type:'image' loads image from file (natural image resized) // 'spikes' creates a synthetic image with k spykes</span>
0023 paramIM.dim=[n m];
0024 paramIM.type=<span class="string">'image'</span>; 
0025 paramIM.k=16;
0026 paramIM.F0=1;
0027 paramIM.fileimage=<span class="string">'eta-carinae_16x16.mat'</span>;
0028 
0029 <span class="comment">% create an image in x</span>
0030 x=gen_astro_object(paramIM); 
0031 <span class="keyword">if</span>(plots), figure, imagesc(x),colorbar, axis image, <span class="keyword">end</span>
0032 
0033         
0034 
0035   
0036 <span class="comment">%% Generate mask for the open frequences</span>
0037 
0038 u=1; <span class="comment">%undersampling ratio u=M/N;</span>
0039 sigmau=1/4;sigmav=1/4;sigmaw=1/4; <span class="comment">%standard deviation of the gaussian profile</span>
0040 T=mask_3D_3(n,m,u,sigmau, sigmav, sigmaw);
0041 Tr=build_redundant_table( T );
0042 
0043 
0044 <span class="comment">%% Add noise only to the non-redundant part of the measurements vector</span>
0045 xhat=1/n*fft2(reshape(x,n,m));
0046 ymeas=Ayz(x(:),xhat(:), xhat(:),T);  <span class="comment">%true measurements without noise</span>
0047 NB=numel(ymeas);
0048 
0049 
0050 <span class="keyword">if</span> noisy
0051     <span class="comment">% Add Gaussian i.i.d. noise</span>
0052     
0053     eB=sqrt(1/NB*sum(abs(ymeas(2:end)).^2));
0054     sigma_noise=.5*10^(-input_snr/20)*eB;
0055     
0056     nBr=(randn(size(ymeas))); 
0057     nBi=(randn(size(ymeas)));
0058     
0059     nBr=sqrt(NB)*sigma_noise*nBr./norm(nBr(:),2);
0060     nBi=sqrt(NB)*sigma_noise*nBi./norm(nBi(:),2);
0061     
0062     nB=nBr + 1i*nBi;
0063      
0064     ymeas=ymeas+nB;
0065     
0066     
0067 <span class="keyword">end</span>
0068 
0069 y=build_redundant_meas(ymeas,T);
0070         
0071 
0072 <span class="comment">%% Parameters of the algorithm</span>
0073 
0074 param.n=n;param.m=m;
0075 param.T=Tr;
0076 param.rel_obj=1e-3;
0077 param.max_iter=1000;
0078 param.verbose=1;
0079         
0080         
0081 <span class="comment">% compute starting estimate for the Lipschitz constant for the backtracking procedure</span>
0082 aux=rand(size(x));aux=aux./sum(aux(:));
0083 auxhat=1/n*fft2(reshape(aux,n,m));
0084 A=@(x)Ayz(x,auxhat(:),auxhat(:),T);
0085 At=@(x)Ayz_t(x,auxhat(:),auxhat(:),T);
0086 
0087 param.L=norm(grad_of( aux, y, A, At ),2)/norm(aux,2);
0088 param.eta=1.1; 
0089 param.max_iter_BT=100;
0090 
0091 max_iter=5; <span class="comment">%max number of time we reinitialize the problem with rand point</span>
0092 of_vals=zeros(max_iter,1);
0093 sol_vals=zeros(max_iter,n,m);
0094 prev_sol=zeros(N,1);
0095 iter=0;
0096 crit=0;
0097 
0098 <span class="comment">%% solve AM problem</span>
0099 t_start=tic;
0100 
0101 <span class="comment">% loop for different initializations of the problem</span>
0102 <span class="keyword">for</span> iter=1:max_iter
0103 
0104     <span class="comment">%initialize x,y,z</span>
0105     param.xinit = rand(size(x(:)));param.xinit=param.xinit(:)./sum(param.xinit(:));
0106     param.yinit=param.xinit;
0107     param.zinit=param.xinit;
0108 
0109 
0110     [sol, of, crit_AM]=solve_am(y,param);
0111 
0112 
0113     of_vals(iter)=of;
0114     sol_vals(iter,:,:)=reshape(sol,n,m);
0115 
0116 <span class="keyword">end</span>
0117 t_end=toc(t_start);
0118 
0119 <span class="comment">% choose the solution corresponding to the minimum objective function</span>
0120 [Max, Max_i]=min(of_vals(:));
0121 sol=sol_vals(Max_i,:,:);sol=sol(:);
0122 <span class="comment">% projection (real)positive orthant</span>
0123 sol=real(sol); sol(sol&lt;0)=0;
0124 
0125 x_AM=reshape(sol,n,m);
0126 
0127 <span class="comment">%% Log</span>
0128 
0129 maxsnr=calculate_snr(x,x_AM);
0130 <span class="comment">% plot the result using the same scale as for the original plot</span>
0131 mx=max(x(:));
0132 <span class="keyword">if</span> (plots&gt;0), figure, imagesc(x_AM,[0,mx]),colorbar, axis image, <span class="keyword">end</span> 
0133 
0134 fprintf([<span class="string">'\n M/N= %e, SNR0 = %e, ellapsed time= %e\n\n'</span>], u, maxsnr, t_end);
0135 
0136 
0137 
0138     
0139 
0140 
0141 
0142 
0143 
0144 
0145 
0146</pre></div>
<hr><address>Generated on Thu 18-Jul-2013 19:25:15 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>